Taylor's swimming sheet: Analysis and improvement of the perturbation series 
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In G.I. Taylor's historic paper on swimming microorganisms, a two dimensional sheet was pro- 
posed as a model for flagellated cells passing traveling waves as a means of locomotion. Using a 
perturbation series, Taylor computed swimming speeds up to fourth order in amplitude. Here we 
systematize that expansion so that it can be carried out formally to arbitrarily high order. The 
resultant series diverges for an order one value of the wave amplitude, but may be transformed into 
series with much improved convergence properties and which yield results comparing favorably to 
those obtained numerically via a boundary integral method for moderate and large values of the 
wave amplitudes. 



I. INTRODUCTION 



In his landmark 1977 paper, Purcell elucidated the unique challenges faced by microorganisms attempting to 
propel themselves in an inertia-less world pQ. In the creeping flow limit, viscous stresses dominate, and thus shape- 
changing motions which are invariant under time reversal cannot produce any net locomotion - the so-called scallop 
theorem. In order to circumvent this limitation many microorganisms are observed to pass waves along short whip-like 
appendages known as flagella, usually transverse planar waves for many flagellated eukaryotic cells, and helical waves 
for prokaryotes PHI]. 

In the first of a series of pioneering papers on the swimming of microorganisms, G.I. Taylor investigated back in 
1951 such motions by considering the self propulsion of a two-dimensional sheet which passes waves of transverse 
displacement [5] . By stipulating that such waves have a small amplitude relative to their wavelength Taylor utilized a 
perturbation expansion to compute the steady swimming speed of the sheet to fourth order in amplitude. Drummond 
later extended Taylor's calculation of the swimming speed of an oscillating sheet to eighth order in amplitude [BJ. A 
concise presentation of the derivation can be found in Steve Childress' textbook [7]. 

Since then many more sophisticated theoretical and computational models have been proposed to study the loco- 
motion of microorganisms which are well documented in several review articles PHI]- Nevertheless, the simplicity of 
the swimming sheet still provides opportunity for insight and analysis into such problems as swimming in viscoelas- 
tic fluids [5] [5] , the synchronization of flagellated cells [TD1 E] , or peristaltic pumping between walls [T^HTB] . The 
swimming sheet has also been utilized to yield theoretical insight into inertial swimming [16j [17] . 

In this paper we show that through some mathematical manipulation the perturbation expansion for an inextensible 
sheet outlined by Taylor (j]TT]) may be performed systematically so that the result may be obtained to arbitrary order 
in amplitude QUI). The resulting series obtained is found to be divergent for order one wave amplitudes. Using 
boundary-integral computations as benchmark results ft IV I, we show however that the series may be transformed 



to obtain an infinite radius of convergence (SjVJ, thus providing an analytical model valid for arbitrarily large wave 
amplitude. The coefficients for both the original and the transformed series are included as supplementary material. 



II. SERIES SOLUTION FOR TAYLOR'S SWIMMING SHEET 



A. Setup 

We consider a two dimensional sheet of amplitude b which passes waves of transverse displacement at speed c — oj/k, 
where ui is the frequency and k is the wavenumber (see Fig. [I]). The material coordinates of such a sheet, denoted by 
s, are given by 

y s = bsin(kx — ut). (1) 

We use the following dimensionless variables for length x* = xk and time t* = tu: (where *'s indicate dimensionless 
quantities). The ratio of the amplitude of the waves to their wavelength is given by e = bk. For convenience we use 
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FIG. 1. Left: Graphical representation of the Taylor's swimming sheet. Right: Example of wave amplitude studied in this 
paper: e = 0.1, 1, 7. 



the wave variable z = x* — t* and therefore write 

y* = esin(z) = ef(z). 



(2) 



The regime we consider here, that of microorganisms, is the creeping flow limit governed by the Stokes equations 
for incompressible Newtonian flows 



V • u* = 0, 
Vp* = V 2 u*, 

where the velocity field u* = {u, v}/c and pressure field p* — p/fiuj. We now drop the *'s for convenience. 
In two dimensions the continuity equation is automatically satisfied by invoking the stream function tp where 



dtp 



dip 



dy ' dx 

The Stokes equations arc then transformed into a biharmonic equation in the stream function 

v 4 v> = o. 



(3) 
(4) 
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The components of velocity of a material point of the sheet are denoted by uq and vq. The conditions to be satisfied 
by the field ip at the surface y = e/(z) are hence 



dip. 



g y \v=ef-U , 

In order to find an analytical solution we seek a regular perturbation expansion in powers of e, 
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where K is the order to which we wish to take our expansion. 

We consider here only the upper-half solution which, by symmetry, is sufficient to yield the swimming velocity. The 
solution to the biharmonic equation which yields bounded velocities in the upper half plane, at 0(e fc ), is given by 



3 = 1 



(Af ] + B^y) sin(jz) + {C) K > + Df'y) cos(jz) 



-jv 



(10) 
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We look to solve this problem in a frame moving with the sheet and hence the terms U^ k 'y allows for the waving 
sheet to move relative to the far field with a velocity equal toU~- J2k=a e k U( k ^e x . 

In order to express the stream function on the boundary we expand ip in powers of e about y = using Taylor 
expansions, and get 



<hp t 
dy 



k-1 



\V=*f 



E^E 
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fit Qn+l^(k—n) 

n\ dy n+ 



\ — lw=o, 



(11) 



and 



dtp . ^ fe f* <9"+ V fe ~ n) 



k=l n=0 



n! dxdy r 



(12) 



Substituting for ip from Eq. (10) and equating with the boundary conditions wc find that for k £ we must 

have 
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At order k, the unknowns, which are the k th coefficients Af\ Bf\ cf\ Df ] and swimming speed U^ k \ are in 
the n = term only. Factoring this off and rearranging we obtain 



(jAf } - SW) sin(jz) + (jCf' - Df>) cos(jz) 



and 



(15) 



E 



jAf cos(jz) - jC (fc) sin(jz) 



with (§W and fiW given by 



fe— 1 k—n 



G (fe) - E E 

n=l j=l 
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(16) 



(17) 
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Provided the solution for the flow field is known for all orders up to k — 1, the left-hand side of Eqs. ( 15 )-( 16 1 is thus 
known, and all the unknowns, determining the fc th order terms, are on the right-hand side. 
The terms and may conveniently be rearranged into a Fourier series of order k as 



k k 

&k) = E k f ] cos (^) + E ^ (fe) sin ^ z )' 



3=0 
k 



3 = 1 
k 



fj(k) = £ f(k) CQs{jz) + J2 Rf sin(jz) 

3=0 j=l 



(19) 
(20) 
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A simple expression of the Fourier coefficients is not easily obtained; however, they are easily (numerically) computed. 

Finally, as we show below, the k th term of the components of velocity at the boundary can be written as a Fourier 
cosine series of order k 

=$>f cos(jz), (21) 

3=0 



and 



,(*) _ «( fe ). 



- - ^PT> cos{jz). (22) 

We can hence write, for all k, the system to solve 

jAf -Bf ] =§f\ (23) 

jCf-Df =afi+kf\ (24) 

jCf^Rf, (25) 

3 Af =$f-ff\ (26) 



for j £ [1, k], or more compactly 



JiAf=hf\ (27) 

The determinant of the coefficient matrix det(Jj) = j 2 and hence invertible Vj 7^ 0. The solutions for each j are 
decoupled, and thus for each k we invert a 4k block diagonal matrix. 
Note that for the mean j = terms we obtain 

tfW = _ tt < fc > _ (28) 

= -f (fc) . (29) 

We thus see that the swimming speed at 0(e fe ) depends only on the mean at that order. We also find that since 

(k) ~(k) 

there is no far-field vertical velocity we require pg = T , which are both known, in order to avoid an ill-posed 
problem. This means that since we do not allow a mean vertical flow in the solution of the stream function (which 

gives Tq = 0) then the vertical boundary conditions must have zero mean, /3q = 0. 

Now we can solve for the swimming speed up to 0(e k ) by solving the above system all k orders sequentially, provided 
we have the Fourier coefficients for the boundary conditions up to 0(e fc ). 



B. Boundary conditions 

Following Taylor [5] , we wish the material of the sheet to be inextensible. In a frame moving at the wave speed the 
shape of the sheet is at rest [SJ [7] , therefore in a frame moving with the sheet the boundary conditions are 

m = -Qcos6 + f, (30) 
v o = —Qsin9, (31) 

where tan 9 = y' s and Q is the material velocity in the moving frame is given by 

1 r 27r 

Q = — V 1 + e2 cos 2 (z)dz. (32) 
2tt Jo 

Expanding in powers of e and integrating we obtain 



2n 



00 

E^ 2 "- ( 33 ) 
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Similarly we expand cos 9 in powers of e to give 

-.-f«-(-iri(*)f-(^) +2 E(^ r )-o») 

oo n 

= E e2 "E^ cos ( 2 ^)- ( 34 ) 

n=0 r=0 

Letting fc — 2n and considering only even values we obtain 

oo fc/2 fe/2 

- 1 - E ek E cos ( 2rz ) E 

fe=0 r=0 p=r 
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= -E efe E cos ( 2r *)E^t- P - ( 35 ) 

fc=2 r=0 p=r 

Now letting j = 2r we find for even j and even fc > 2 



fc 

( fc ) - ^,( fc ) 
o 



E a f cos(jz), (36) 



where we have defined 

.,0) 



fe/2 

- E ^ 9 f-p' (37) 

P=j/2 



while = for odd fc and af^ = for odd j. 
We then know that 

v Q = -y' s (z)Q cos 6, (38) 

and hence we find for odd j and odd k > 3 



( fe ) - \^ «« 



£$ fc) cos(j*)> (39) 

= at 1] + (40) 

(fc-i) (fc-i) 

ti k) = j ~ l 2 j+1 , 3<i<fc-2, (41) 

= (42) 

and for fc = 1 fi^ = —1. In contrast, Vq = for even fc and = for even j. We see that the vertical component 
of the boundary velocity has no mean component at any order in e which, as we saw in the previous section, is required 
given the form of the solution. With these coefficients we can now solve a linear system at each order to obtain L/w 
to arbitrary order. 

In practice the number of terms obtainable will be limited by numerical technique. To obtain the first one thousand 
terms of the series used in the analysis in the following sections, the system of equations was solved using the 
C programming language with GNU MP, the GNU multiple precision arithmetic library [TH], using 300 digits of 
accuracy. 

III. ANALYSIS AND IMPROVEMENT OF THE PERTURBATION SERIES 

In the previous sections we presented methodology to obtain the solution to the swimming speed U in the form of 
a perturbation series 

K 



C/(e)-^[/ (fe) e fc . (43) 



fc=i 



G 




It remains of course to be seen whether the series will converge to U for arbitrary e. We analyze here the convergence 
properties of the series, and methods to improve upon that convergence. 

In Fig. [2] we plot the coefficients of the series against k. On Fig. [2]} are plotted the first 100 terms, and on 
Fig. [2]d the logarithm of the absolute value of the nonzero terms up to k = 1000. We see that the coefficients have an 
exponentially increasing amplitude while alternating in sign, U (k) > for k = An - 2 and Z7 (fc) < for k = An where 
neN. We also note that due to the e — > — e symmetry of the geometry in the problem, all odd powers in the series 
are zero. It is therefore useful to recast the series as follows 

K K/2 K/2-1 

U = J2u^e k = J2u^e 2k = 6 £ c k S k , (44) 

k=l fc=l k=0 

where = [/( 2fc + 2 ) and S = e 2 . The coefficients for k — 1 to 500 have been reproduced in the included supplemen- 
tary material of the manuscript. 

The sign of c k alternates in a regular manner which indicates that the nearest singularity lies on the negative real 
axis and since only positive values of S have any meaning, there is no physical significance to the singularity; it does 
of course govern the radius of convergence of the series |19j . 



A. Series convergence 

The radius of convergence, do, of the power series 

/(,5)~5> fc 5 fc , (45) 

k 

may be simply found by using the ratio test 

60 = lim — ■ (46) 

k— f 00 C/c 

In order to find this value we must extrapolate due to the finite number of terms. In order to aid this process Domb 
and Sykes noted it is helpful to plot Ck/ck^i against 1/k [2D]. The reason is that if the singular function /, has a 
dominant factor 

(6 -6y for 7^0,1,2,..., (47) 
(6 - £) 7 m(5 - 5) for 7 = 0,1, 2,..., (48) 
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FIG. 3. Domb-Sykes plot of the coefficients Ck, from the series in Eq. ( 44 1 , shows convergence to l/<5o ~ —1.093 



then the coefficients behave like 



Cfc 



1 + 7 



(49) 



for large k [23 [ST]. The result in Eq. (49) indicates that the intercept l/k = in a Domb-Sykes plot gives the 
reciprocal of the radius of convergence while the slope approaching the intercept gives 7. In Fig. [3] we show the 
Domb-Sykes plot of the series Cfc. The plot indicates that the nearest singularity is at Sq « —0.914912217581184, and 
that 7 = — 1 corresponding to a first order pole. 



B. Euler transformation 

One approach to improve convergence of the series is to factor out the first-order pole characterized above, and 
characterize the singularities of the new series. However we find that that series is no more tractable due to the 
presence of an apparent branch cut in the complex plane close to 5 = — 1 . 

Alternatively, the original non-physical singularity Sq may be mapped to infinity using a Euler transformation and 
introducing a new small variable 



The power series for / is then recast as 



S = 



S-6 



(50) 



(51) 



The coefficients du for k = 1 to 500 have been reproduced in the included supplementary material of the manuscript. 
Their values for k > 50 are shown in Fig. [4^i and we can see that they decay in magnitude for large k. 

In order to find the radius of convergence of the new series, dk, we again turn to the Domb-Sykes plot, which is 
show in Fig. [4J5. We see that it appears dk/dk-i — > 1 as fc _1 — > and since S/(S — <5o) — > 1 when 5 — > 00, we have 
now have an infinite radius of convergence in the original variable S. Note that the vastly improved convergence does 
not necessarily mean the series will provide a good approximation beyond <5 |21j ; however we will see in the results 
section that it actually provides an excellent fit to the numerical results. 




C. Pade approximants 



A popular scheme to improve the convergence properties of series is to recast the series as a rational polynomial 

M+N X k 

m „ y: ^ ~ = P M t (52) 

fc=0 



where M + N < K/2 — 1. If we multiply both sides by the denominator Yl bkS k for the terms of order S k where 
k = M + 1 : M + N we obtain a square matrix to invert for bi, and one takes bo = 1 with no loss of generality 
|22j . One can then solve for afc. 

We apply this method to our swimming sheet, and plot the zeros of different Pade denominators with M = N in 
Fig. [5] It is evident that the pole we identified earlier at So is well reproduced here. The interesting feature beyond 
this is the fact that the remaining zeros do not exhibit any consistency, which indicates a branch cut in the complex 
plane. 

D. Shanks transformation 

A scheme to improve the rate of convergence of a sequence of partial sums 

n 

S n = Y,CkS k , (53) 

fe=0 

for n — to < K/2 — 1, is to assume they are in a geometric progression 

S n = A + BC n . (54) 
Solving for A by nonlinear extrapolation of three sums yields 

4 _ q (ffn+i — S n ) (S n — S n -i) , . 

The A n 's for n = 1 to N — 1, can then be considered a series of partial sums and the Shanks transformation may be 
thereby repeated (N — l)/2 times [2Tj . 
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FIG. 5. Zeros in the complex plane of the denominators of various Pade approximants for M = N = 10, 50, 100, 200 and 249. 



IV. BOUNDARY INTEGRAL FORMULATION 



In order to provide benchmark results for the analysis of the perturbation series and its various transformations, 
we use the boundary integral method to obtain what we will consider to be an exact solution of the swimming speed 
for waves of arbitrarily large amplitude. 

We briefly summarize the principle of the method here. The Lorentz reciprocal theorem states that two solutions 
to the Stokes equations, (u, <x) and (u, <x) are related by 

/ (u • a) ■ n dS = [ (u-a-)-ndS, (56) 
Js Js 

within a volume V bounded by the surface S whose unit normal n is taken pointing into the fluid. The velocity and 
stress fields, u(x) and <r(x), are taken to be fundamental solutions for two-dimensional Stokes flow due to a point 
force at xo, 

u(x) = i-G(x) • f (x ), (57) 

or(x) = ^T(x) • f (x ), (58) 

where x = x — Xo and the two dimensional Stokeslet G, and stresslet T are given by 

G = -Iln(|x|) + -_ (59) 
l x l 

T = -4^J- (60) 
x 



Taking the singular point x to be on the boundary S one obtains from Eq. ( 56 ) a boundary integral solution to 
two-dimensional Stokes equations for the velocity 

u(x ) = A. \ (u(x) • T(x) • n(x) - f(x) • G(x)) dS(x), (61) 
2tt Js 
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where f = cr ■ n. 

We wish to capture the swimming speed of an infinite sheet therefore the domain of integration is an entire half 
plane of fluid bounded by the sheet. In order to avoid performing an integration over the entire bound it is convenient 
to use an array of periodically placed Stokeslets and stresslets, given by 



OO 

E 

n=— oo 
oo 

71 — — OO 



-Hn(|x„|) 



v n- A -n- A -n 



(62) 
(63) 



where x„ = {£o + 27m, z/ }, so that we may then instead integrate G p and T p over a single period [23]. The 
periodic Stokeslet and stresslet may be conveniently expressed in closed form [231 HI] > through the use of the following 
summation formula 



.4 



OO 
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n=— oo 



ln(|x n |) = i In [2cosh(y ) - 2cos(x )] 



(64) 



and its derivatives, as follows 
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yy 



-A - d y A + 1, 
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and 



rpp 

xxy 

Txyy = tydxyA, 
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-2d x (2A + yd v A), 
-2d y (yd y A), 



~2{dyA - ydyyA) 



(68) 
(69) 
(70) 
(71) 



The remaining elements follow from a permutation of the indices of the Stokeslet and stresslet which leaves the right 
hand side unchanged fM\ . 

The flow is quiescent at infinity and periodic on 2n and therefore the domain of integration S reduces to the surface 
of the sheet over one period. To facilitate integration the continuous boundary is discretized into N straight line 
elements S n and we assume that f is a linear function over each particular interval, f — > f n (see Ref. We 
decompose the boundary velocity into surface deformations and rigid body motion u — > u„ + U, where u„ is a linear 
function over each interval and U = —Ue x . Then xo is taken at the center of each of the the N segments S n ,where 
the velocity is known, xo — > x m . The G p and T p are regularized by subtracting off the Stokeslet and stresslet from 
their periodic counterparts. The two-dimensional Stokeslet and stresslet are then integrated analytically and added 
back. 



We thereby obtain from Eq. (61 1 a linear system for f„ and U, given by 



u(x m ) + U = 



1 N 

— Y 

2?r ^ 

n=l 



f f„ • (G p - G) dS n - f f„ • GdS n 



K + U) • (TP - T) • n n dS n 



(u„ + U) • T ■ n n dS n 



(72) 



We then obtain U by specifying that the sheet is force free 

N 



E 



f n dS n 



Sn 



(73) 



The numerical procedure was validated by reproducing Pozrikidis' results for shear flow over sinusoidal surface [23] . 
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FIG. 6. Swimming speed, U, against wave amplitude, e, for the unaltered series, Eq. (43 1, with K = 4 (dashed-dot), K = 20 
(dashed), K — 1000 (solid). The series diverges for e « 0.9565. Red squares indicate data points from the boundary integral 
method. 



COMPARISON BETWEEN SERIES SOLUTION AND COMPUTATIONS 



Series solution 



We first show the convergence of the unaltered series expansion, Eq. (43 1, in Fig. [6] where we display the swimming 
speed of the sheet, U, as a function of its amplitude, e. The red squares indicate numerical points computed with the 
boundary integral method. We plot the results for Taylor's original fourth order expansion (dashed-dot line) which is 
reasonably accurate up to e ps 0.4. The series with iiT = 20 is shown in dashed line. As we add terms, we get that the 
series with K = 1000 (solid line) fails to converge beyond the singularity at e = y/—5a ~ 0.95651, as expected from 
the analysis in ^ III 



Euler transformation 



The presence of the singularity on the negative real axis for the series Cfc led naturally to an Euler transformation 
to map the singularity to infinity which, as detailed in §111 B[ yields a series with an infinite radius of convergence in 
8, and thus in e. In Fig. [7] we plot the results of the Euler-transformed series, Eq. (51), for the swimming speed, U, 
against the wave amplitude, e. The results are markedly improved over the original unaltered series. With K = 4 
we obtain results which are accurate for up to e ~ 1.3, already higher than for Taylor's fourth order formula. With 
K = 20 terms, U(e) is found to be accurate up to e ~ 2, and when using K — 100 terms we obtain results which are 
accurate for e > 7. With all K = 500 terms the series is accurate up to e ~ 15 with a relative error of 1% (the series 
is however convergent for all values of e). 



C. Pade approximants and Shanks transformation 



Pade approximants provide a convenient (yet brute-force) way to drastically improve the performance of the series 
without the need to investigate the analytic structure of the underlying function. We find that using only a few 
terms provides very good results, as we show in Fig. ^1 For K = 4 we obtain P| (dashed) which is accurate past 
the singularity, while for K = 22 we obtain P^q (solid) which is accurate up to e w 4, and shows an error which is 
reasonably small for larger amplitudes. Unfortunately the coefficient matrix which must be inverted to obtain the bk 
coefficients of the Pade approximants becomes increasingly ill-conditioned as more terms of the series are added and 
we see diminishing returns from the Pade approximants of higher order expansions; for example, is only accurate 
up to £ K 5. 

Similarly, repeated Shanks transformations of the first few partial sums results in a marked improvement of the 
convergence of the series. We see in Fig. [8] that the (repeated) Shanks transformation of partial sums up to £2 (dotted 
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FIG. 7. Swimming speed, U, against wave amplitude, e, for the Euler series, Eq. (51 1, with K = 4 (dashed-dot), K — 20 
(dashed), K — 100 (solid). Red squares indicate data points from the boundary integral method. 
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FIG. 8. Swimming speed, U, against amplitude, e, for (repeated) Shanks transformations of partial sums up to S2 (dotted) and 
56 (dashed-dot) and for the Pade approximants P| (dashed), (solid). Red squares indicate data points from the boundary 
integral method. 

line) yields results nearly identical to the P| approximant, while for terms up to 5*6 (dashed-dot) we see reasonable 
accuracy up to e ~ 2 in agreement with the results from Ref. [5J. We find however that the addition of any further 
terms in the sequence leads to a pronounced decrease in the convergence properties of the sum. 



VI. CONCLUDING REMARKS 



Despite its simplicity, Taylor's swimming sheet model is still used to provide physical insight into many interesting 
natural phenomena. In this paper, we demonstrated that by systematizing the perturbation expansion outlined by 
Taylor in the wave amplitude, e, the solution for the swimming speed can be obtained in a straightforward fashion to 
arbitrarily high order. The series unfortunately diverges for e rs 0.9565 due to a nonphysical first-order pole located 
in the negative real axis. In order to increase the convergence of the series, the singularity can be mapped to infinity 
via an Euler transformation. The recast series then has an infinite radius of convergence and produces spectacularly 
accurate results for very large amplitudes (albeit requiring a good number of terms). An alternative is to reformulate 
the series using Pade approximants or repeated Shanks transformations, which give reasonable accuracy for moderate 
amplitudes with fewer terms, but can become problematic for very large amplitudes. 
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